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The usual explicit finite-difference method of solving partial differential equations is limited in 
stability because it approximates the exact amplification factor by power-series. By adapting the 
same exponential-splitting method of deriving symplectic integrators, explicit symplectic finite- 
difference methods produce Saul'yev-type schemes which approximate the exact amplification factor 
by rational-functions. As with conventional symplectic integrators, these symplectic finite-difference 
algorithms preserve important qualitative features of the exact solution. Thus the symplectic diffus- 
ing algorithm is unconditionally stable and the symplectic advection algorithm is unitary. There is a 
one-to-one correspondence between symplectic integrators and symplectic finite-difference methods, 
including the key idea that one can systematically improve an algorithm by matching its modified 
Hamiltonian more closely to the original Hamiltonian. Consequently, the entire arsenal of symplec- 
tic integrators can be used to produce arbitrary high order time-marching algorithms for solving the 
diffusion and the advection equation. 



I. INTRODUCTION 

The ID diffusion equation 

can be solved numerically by applying the forward-time and central-difference approximations to yield the explicit 
algorithm 

u'j = Uj + r[(u j+ i - Uj) - (Uj - Uj-i)], (1.2) 
where Xj — jAx, Uj = u(xj,t), u'j — u(xj,t + At) and 

AtD , . 

r = — T . (1.3) 
Ax 

Under this (Euler) algorithm, each Fourier component Uk = e lkx with wave number k is amplified by a factor of 

g = 1 - 4rsin 2 (fcAr/2), (1.4) 
restricting stability (|<7| < 1) to the Courant-Friedrichs-Lewy-i (CFL) limit, 

r<\- (1-5) 



Since explicit finite-difference methods approximate the exact amplification factor by power-series such as (|1.4I) . it 
seems inevitable that they will eventually blow-up and be limited in stability. However, Saul'yeviLul showed in the 50's 
that, by simply replacing in (11.21) . either 

( Uj - Uj-x) (u^ - Uj-_i) or (u j+1 - Uj) (u' j+1 - u'j) (1.6) 

one would have unconditionally stable algorithms: 

u'j = Ps Wj-i + 7s Uj + [5 S Uj + i, (1.7) 

or 

u 'j = Ps Uj-i + 7s Uj + I3 S u' j+1 , (1.8) 
where 7s and (3s are Saul'yev's coefficients given by 

7s = -7-T~ and Ps = 7~ — ■ (!- 9 ) 
1 + r 1 + r 
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Algorithm (|1.7j) is explicit if it is evaluated in ascending order in j from left to right and if the left-most u\ is a 
boundary value fixed in time. Similarly, algorithm (|1.8I) is explicit if it is evaluated in descending order in j from right 
to left and if the right-most un is a boundary value fixed in time. Saul'yev also realized that both algorithms have 
large errors (including phase errors due to their asymmetric forms), but if they are applied alternately, the error would 
be greatly reduced after such a pair-wise application. This then gives rise to alternating direction explicit algorithms 
advocated by Larkin^ and generalized to alternating group explicit algorithms by Evans^. 

There are four unanswered questions about Saul'yev asymmetric algorithms: 1) While it is easy to show that 
algorithm (|1.7[) and (|1.8I) are unconditionally stable, there is no deeper understanding of this stability. 2) The 
algorithms are not explicit in the case of periodic boundary. What would be the algorithm if there are no fixed 
boundary values? 3) The alternating application of (|1.7p and (11.81) greatly reduces the resulting error. How can one 
characterize this improvement precisely? 4) How can Saul'yev-type algorithms be generalized to higher orders? 

This work presents a new way of deriving finite-difference schemes based on exponential-splittings rather than Taylor 
expansions. Exponential-splitting is the basis for developing symplectic integrators^—, the hallmark of structure- 
preserving algorithms. The finite-difference method presented here is properly "symplectic" in the original sense that 
it has certain "intertwining" quality, resembling Hamilton's equations. It is also symplectic in the wider sense of 
structure-preserving, in that there is a Hamiltonian-like quantity that the algorithms seek to preserve. 

As will be shown, there is a one-to-one correspondence between symplectic finite-difference methods and symplec- 
tic integrators. It is therefore useful to summarize some basic results of symplectic integrators for later reference. 
Symplectic integrators are based on approximating e e ( A+B ) to any order in e via a single product decomposition 

ee( A+B) = -Q e a l eA e 6 l eB ; (1 1Q) 

i 

where A and B are non-commuting operators (or matrices). Usually, A+B=H is the Hamiltonian operator and e eH 
is the evolution operator that evolves the system forward for time e. The key idea is to preserve the exponential 
character of the evolution operator. The two first-order, Trotter— approximations are 

T 1A (e) = e e V B , T 1B ( e )=e eB e eA , (1.11) 

and the two second-order Strang-^ approximations are 

T 2A (e) = r 1A (e/2)T 1B (e/2)=e^ A e £B e^ A , 

T 2B (e) - T 1B (e/2)T 1A (e/2)=e^ B e eA ei £B . (1.12) 

The approximation 

T 2C (e) = ^[r 1A ( e )+T 1B (e)] (1.13) 

is also second-order, but since it is no longer a single product of exponentials, it is no longer symplectic. In most 
cases, it is inferior to T 2A and T 2 b because the time steps used in evaluating Ti A and Tib are twice as large as those 
used in T 2A and T 2B . 

Let T 2 denotes either T 2A or T 2B . T 2 must be second order because for a left-right symmetric product as above, it 
must obey 

and therefore must be of the form, 



T 2 (-e)T 2 (e) - 1, (1.14) 



T 2 (e) = '•• (1.15) 

with only odd powers of e in the exponent. (Since there is no way for the operators in (|1.14[) to cancel if there are any 
even power terms in e.) In (| 1 . 15|) . E„ denote higher order commutators of A and B. The algorithm corresponding to 
T 2 then yields exact trajectories of the second-order modified Hamiltonian 

H 2 (e) =H + e 2 E 3 + e 4 E 5 + ---. (1.16) 

A standard way of improving the efficiency of symplectic integrators is to generate a 2n*' t -order algorithm via a 
product of second-order algorithm o 7 ' 9 ' 10 , via 

N 

T2n{e) = Y[T 2 (a l e). (1.17) 
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Since the error structure of 12(e) is given by (|1.15[) . to preserve the original Hamiltonian, one must choose a, to 
perserve the first power of e, 



N 



J2 a > = L 



(1.18) 



To obtain a fourth-order algorithm, one must eliminate the error term proportinal to e 3 by requiring, 



N 



(1.19) 



For a sixth-order algorithm, one must require the above and 

N 



(1.20) 



and so on. While proofs of these assertions in terms of operators are not difficult, we will not need them. Symplectic 
finite-difference methods use a much simpler version of these ideas. Instead of dealing with the evolution operator 
e eH , the finite difference method has a proxy, the amplification factor, which is just a function. Order-conditions such 
as (|1.19p and (|1.20j) will then be obvious. Other results will be cited as needed, but these basic findings are sufficient 
to answer the four questions about Saul'yev's schemes. For the next two sections we will give a detailed derivation of 
the symplectic diffusion and advection algorithms, followed by a discussion of the diffusion-advection equation and a 
concluding summary. 



II. SYMPLECTIC DIFFUSION ALGORITHM 



Consider solving the diffusion equation (jl.ip with periodic boundary condition ujv+i = u\ in the semi-discretized 
form, 



Regarding Uj as a vector, this is 



with 



duj D 

~dt = Ax* ( Uj+1 ~ 2U] + Uj - 1 '- 



du 

-dt =AU ' 



D 

A? 



(-2 1 
1 -2 1 



V 1 



1 -2 1 
1 -2/ 



and exact solution 



u{t + At) =e AtA u{t). 

The Euler algorithm corresponds to expanding out the exponential to first order in At 

u(t + At) = (1 + AtA)u{t), 

resulting in a power-series amplification factor (jl.4p . with limited stability. 

If the exponential in (|2.4I) can be solved exactly, the amplification factor would be 



(2.1) 



(2.2) 



(2.3) 



(2.4) 
(2.5) 



Qex — ^ 



(2.6) 
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where 



r4sin 2 (6>/2) and 9 = kAx. 



(2.7) 



The amplification exponent h ex here plays the role of a time parameter r times the original "Hamiltonian" ho = 
4sin 2 (6*/2). The resulting algorithm will then be unconditionally stable for all r > 0. In the limit of Ax — > 0, each 



/c-Fourier components will be damped by g e 



-AtDk* 



which is the exact solution to 



To preserve this important feature of the exact solution, one must seek alternative ways of approximating of e AtA 
without doing any Taylor expansion. The structure of A immediately suggests that it should decompose as 



/V 



(2.8) 



where each Aj has only a single, non- vanishing 2x2 matrix along the diagonal connecting the j and the j + 1 elements: 



A,- = 



D 
A? 



/. 



V 



-1 1 
1 -1 



/ 



and A^v = 



The exponential of each Aj can now be evaluated exactly: 



1 



.AtAi 



a (3 
(3 a 



D 
A? 



/-l 



V 1 



1 \ 



-1/ 



(2.9) 



(2.10) 



where 



a = ~(l+7), P = \0- 



7), and 7 



(2.11) 



2" 1,1 ' 2 
Each e AtAj updates only Uj and Uj+i as 

u'j = au.j + (3uj+i 

u'j + i = /3uj+auj+i. (2.12) 

The eigenvalues of this updating matrix are a ± /3 = 1, 7, with det = 7. This means that the updating is dissipative 
for r > and unstable for r < 0. Since a and (3 are given in terms of 7, the resulting algorithm depends only on a 
single parameter 7. 

One can now decompose exp(ZitA) to first order in At (apply (|1.11[> repeatedly) via either 

T 1A (At) = e AtAN ■ ■ ■ e AtA *e AtAl , (2.13) 



T 1B (At) = e ntAl ■ ■ ■ e ^»-i e ^». (2.14) 

These algorithms update the grid points sequentially, two by two at a time according to (|2. 12|) . but each grid point 
is updated twice, in an intertwining manner. This is crucial for dealing with the periodic boundary condition. Let 
Uj denotes the first time when Uj is updated and u'j the second (and final) time it is updated. One then has for 
algorithm 1A: 

(2.15) 



* 


= au\ - 


Vfiu 2 


* 

"2 


= /3ui - 


V au 2 




= au* 2 - 




* 

u 3 


= /3«2 - 


- au 3 



U 3 



aiij_|_i. 



'N 



au* N + f3u\ 
(3u* N + au[. 



(2.16) 
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Since a + ft = 1, summing up both sides from (|2. 15|) to (|2.16[) gives, 

N 

4 = I>- ( 2 - 17 ) 



N N 

i=i j'=i 



The algorithm is therefore norm-conserving. The same is true of algorithm IB below. For 2 < j < N one has 



au* + fiuj + i 

a(f3u*_ 1 + otUj) + Puj+i 

f3(u' J _ 1 — f3uj) + a 2 Uj + (3u j+ i 

/3u' j _ 1 +'yuj+l3u j+ i (2.18) 



and for j = 2, N, 



"2 

u' N 

Finally when the snake bits its tail, one has 



Similarly, IB is given by 



'■jv 



j3ul + "fu 2 + (3u 3 

Pu' N _ x +7U N + Pul. (2.19) 



ryft 

4 + 7«i + — «a- (2-20) 
a 



PiiN + aux (2.21) 



u'j = Puj-x + -yuj + j9tt^ +1 (2.22) 
it' 2 = /9uJ + 7712 + /?M3 

u'i = ——u jy + 7U1 + — Uo. (2.23) 
a a 

Algorithms 1A and IB are essentially given by (|2.18j) and (|2.22|) respectively, except for three values of u[, u' 2 
and u' N . The forms of (|2.18p and (|2.22l) reproduce Saul'yev's schemes (|1.7p and (|1.8p . but with different coefficients. 
Saul'yev's coefficient 7s is a rational approximation to the 7 = e~ 2r here. Note that his /3s = r/(l + r) is also given 
by /3s = (1 — 7s)/2- In contrast to Saul'yev's algorithm, which cannont be started for periodic boundary condition, 
algorithm 1A and IB are truly explicit because they are fundamentally given by the sequential updating of (|2.12l) . 
Each algorithm can get started by first updating u\ to u*, then updating it again at the end to u[. 

By virtue of (|1.12p one can now immediately generate a second-order time-marching algorithm via the symmetric 
product, 

T 2 (At) = T 1B (y)T 1A (y) 

_ e \AtK^ e \AtKo. _ _ , e \AtK N e \AtK N , , _ & \AtKi & \AtKi_ (2 24) 

If the boundary effects of u[, u' 2 and u' N are ignored (for now) and 1A and IB are considered as given by (|2.18|) 
and (|2.22p . then the alternative product T\x{At j "£)Txa{At / '2) yields the same second-order algorithm. In this case, 
algorithms 1A and IB have amplification factors 



7 + /3e 11 
1 - /3e- 
7 + /3c" 



■9ib = 1 _ p c +ie . (2-25) 



with opposite phase errors, and the second-order algorithm has 

'/4A /Mi\ 

7 2 + /3 2 + 2/37cos6» 



'/2 = !/i.L'< \ — J giA I -y I 



l + /3 2 -2/3cos6> 



(2.26) 



l-(4/37/5 2 )sin 2 fl/2 = ^ 2 
l + (4/3/5 2 )sin 2 6>/2 
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with no phase error and where 

5=1(1 + 7), £=1(1-7) and 7 = 7 (r/2). (2.28) 

Since both algorithms 1A and IB have phase errors, only the second-order algorithm is qualitatively similar to the 
exact solution. Eq. (|2.27[) makes it clear that this algorithm is unconditionally stable since < 7 < 1. Algorithms 1A 
and IB are also unconditionally stable since |<7ia.ib| = \fgi with 7 — > 7. Note that this also proves the unconditional 
stability of Saul'yev's algorithms. His coefficient 7s can turn negative, but only approaches -1 as r -> 00. Conventional 
explicit methods, like that of the Euler algorithm, are limited in stability because they have power-series amplification 
factors. By contrast, symplectic finite-difference methods are unconditionally stable because they produce SauPyev- 
type schemes with rational-junction amplification factors. This type of stability is usually associated only with implicit 
methods. 




FIG. 1: Comparing the amplification factor of various diffusion algorithms at r = 0.3 and r — 2. Euler is the first-order explicit 
algorithm. CN is the second-order implicit Crank-Nicolson algorithm. D2 is the second-order symplectic algorithm using the 
original coefficients (|2.28[) and D2S is the same alogrithm but uses Saul'pev's coefficient (|2.33[) . "Exact" is g EX of (|2.6[) . 



Since the order of matrices defining T2(At) in (|2.24|) is left-right symmetric, one has the same situation as in the 
symplectic integrators case of (|1.14|) . implying that g2(—r)g2(r) = 1. This means that 172 (— r) > 1 for all k 7^ modes 
and the algorithm is unconditionally unstable for negative time steps. Moreover, this also means that h2(—r) = —fi2(r) 
and /i2(r) is an odd function of r. This is a simpler, functional version of (jl.151) . Expanding /12 of (|2.27p in powers of 
8 gives, 

1+7 6(l + 7) 3 v ; 

Each coefficient must be an odd function of r. This is satisfied only if 

7(-r)7(r) = 1. (2.30) 
Thus every function 7(7") satisfying (|2.30[) with | 7 (r)| < 1 defines an unconditionally stable algorithm for solving the 



7 



diffusion equation. For the original choice of 7(7*) = e r , one finds 

Comparing this to the expansion of the exact amplification exponent, 

h ex = r9 2 - —9 4 + —9 6 + ■■■, (2.32) 
12 360 v ' 

one sees that the original choice does not reproduce leading term exactly except when r << 1. To improve this, let's 
take 7(r) to be an arbitrary function of r but with a and /? still defined by (12 .28[) . The first term in h ex can now be 
matched exactly by requiring 

1 + 7(7") 2 1 + r/2 

which is precisely Saul'yev's original coefficient. With this choice for 7(7*), (|2.27l) reads 

_ l-2r(l-r/2)sin 2 (fl/2) 
92 l + 2r(l+r/2)sin 2 (#/2)' 1 ' 

with exponent 

3 35 

which is now correct to third-order in 9. Comparing this and (|2.31[) to h ex , one sees that all the error terms of h 2 are 
odd powers of r higher than the first. As a matter of fact, by resumming terms proportional to r, we can make this 
error structure in exact conformity with (|f .151) . 

h a = rh + r 3 E 3 (9) + r 5 E 5 (9) + ■■■, (2.36) 

where ho = 4 sin 2 (9/2). This is the same error structure exploited by symplectic integrators to produce higher order 
algorithms. 

Saul'yev's algorithm is close in reproducing the amplification factor of the implicit Crank-Nicolson (CN) scheme 
(which is without the ±r/2 terms in (|2.34jl ). The CN scheme has the advantage that its exponent is 

^ = ^-^ + (4 + ^ + ... (2.37) 

which is correct to fifth-order in 9. In FigJTJ we compare the amplification factor of various algorithms at two values 
of r. D2 and D2S are second-order symplectic diffusion algorithms described above using the original coefficient (|2.28[) 
and Saul'pev's coefficient (|2.33[) . respectively. For r = 0.3, both D2 and D2S track g ex closely over the entire range 
k values. Both Euler and CN tend to over-damp higher Fourier modes. At r = 2.0, while both Euler and CN turn 
negative at large k, D2 and D2S remain positive, like that of g ex . At large r, D2S is clearly better than D2 at small 
k. 

To generalize Saul'yev schemes to periodic boundary condition, one simply replaces in the above algorithms, 7 — > 75. 
In Fig[2]we show the working of algorthms 1A, IB and 2 using both sets of cofficients. With the original coefficients, 
the phase errors are much smaller, but the under-damp error is much larger. For Saul'pev's coefficient, the phase 
errors are much greater, but the under-damp error is smaller. Since the phase error is automatically eliminated by 
going to second-order, D2S has an advantage over D2. 

If one were to construct a 2n" i -order algorithm out of a product of second-order algorithms 

JV 

T 2n (At) = l[T 2 (a t At), (2.38) 

i=i 

then the corresponding h 2n is given by 

N 

h 2n {r) = ^ h 2 (a.ir) 

i=l 

JV JV JV 

= rhoY,^ + r 3 E 3 (9)J2^ + r 5 E 5 (9)J2^ + ---, (2-39) 



8 




FIG. 2: The diffusion of a Gaussian profile after t — 1 on a grid of 120 points spanning the interval [-6,6] with D = 1/2 and 
At = 0.1, corresponding to r = 5. This large value is choosen to exaggerate various errors. The open and filled circles are 
algorithms 1A, IB and 2 using the standard coefficients. The open and filled diamonds are the same three algorithms using 
Saul'pev's coefficients. The solid black line is the exact solution in the continuum limit. 



and the order conditions (jl . 18[) - (II .20[) are easily understood. Unfortunately, for the diffusion algorithm, T 2 (aiAt) is 
unstable for any negative cij and no negative coefficient ai can be allowed. In this case, order-conditions such as (|1.19[) 
and (|1.20p cannot be satisfied and no higher-order composition algorithms of the form (|2.38[) is possible. (However, 
these algorithms will be useful for solving the advection equation in the next section.) 

To overcome this impass, one must go beyond the single product approximation of (|2.38p . and consider a multi- 
product expansion of the form 

N k 

T 2n (At) = J2ckY[T 2 (a kl At). (2.40) 

k i=l 

If a,ki were to remain positive, then Shangi^ has shown that any single product in (|2.40p can at most be second order 
and that some Cfe coefficients must be negative. More recently this author realize tha t 17 ' 18 , due to the error structure 
(|2.36p . the first-order term rho is automatically preserved by monomial products of the form T^At/k) with exponent 

h ( 2 k} = rh + k- 2 r 3 E 3 (9) + k~ 4 r 5 E 5 {8) + ■■■. (2.41) 

The arbitrariness in Nk and can be eliminated by taking Nk — k and a k i — l/k. This then produces a much 
simpler Multi-Product Expansion 17 (MPE) 

T 2n (At) =^c fc T 2 fe (^/fc) (2.42) 

k 

for any sequence of n whole numbers {k} with analytically known coefficients Ck- For the harmonic sequence of 
k = 1, 2, 3, • • •, the first few higher order algorithms are: 



T^t) = - l -T 2 {At) + \Ti{^ 



(2.43) 
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T , m . __L^ + ^ (*) _ ™ rf (-) + ^ (-) . 

For the diffusion equation, the use of the fourth-order extrapolation (|2.43l) has been previously suggested by 
Schatzman 19 . However, these high order methods are uselss for conventional explicit schemes, since they are un- 
stable at large time steps. It is only with the use of unconditionally stable algorithms here that the power of these 
high order schemes can be unleashed. These MPE algorithms do not preserve the positivity of the initial profile. This 
is in keeping with the general observation that there can't be any finite-difference scheme for solving the diffusion 
equation that preserves positivity beyond the second-order—. More recently, Zillich, Mayrhofer and Chin22 have 
shown that Path-Integral Monte Carlo simulations, where positivity is of the utmost importance, can be successfully 
carried out using these expansions, demonstrating that the violation of positivity is small and controllable. These 
MPE algorithms are also not symplectic. However, as argued by Blanes, Casas and Ros^S., they are symplectic to 
order In + 3. Thus at sufficiently high orders, they are indistinguishable from truly symplectic algorithms up to 
machi ne p recision. In the context of classical dynamics, these MPE algorithms have been tested up to the 100 th order 
in Refll8l 




FIG. 3: The convergence of (|a;|) after a Gaussian profile has been diffused for t = 1 on a grid of 120 points spanning the 
interval [-6,6] with D = 1/2. The range of At, from At — 1/400 to At — 1/9, corresponds to a range of r = 0.125 to r = 5.555. 
Lines are fitted power laws At n verifying the order of the algorithm. "12Euler" labels results of running the first order Euler 
algorithm 12 times at time step At/12. "6D2S" are results from running the second-order algorithm D2S six times at time 
step At/ 6. Ti and Tq are fourth and sixth-order algoirthms which require 3 and 6 runs of D2S respectively. 



In FiglSl we compare and verify the order of convergence of T4 and T§ with D2S as T2 by computing the expectation 
value 

(M) = %J^* (2.46) 
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after evolving Uj for t = 1 as a function of At. The absolute value is used because the Euler algorithm would exactly 
preserve (x) even when it is unstable. For computing the Euler algorithm is extremely linear within its tiny 

range of stability. Since it tends to over-damp, its evolving profile is flatter and (|x|) converges from above. Symplectic 
diffusion algorithms under-damp, and their results for (|a;|) converge from below. All results can be well fitted with 
power laws of the form a + bAt n with n — 1,2,4 and 6, verifying the order of the algorithms. To show that these 
higher order algorithms are more efficient than running low order algorithms at reduced step sizes, we also plotted 
results of running the Euler algorithm 12 times at time step At/12 and algorithm D2S six times at At/6. 



III. SYMPLECTIC ADVECTION ALGORITHM 



For the advection equation 



its usual semi-discrete form is 



with discretization matrix 



du 



du 
dx' 



diij i 
dt 2Ax y 1+ 3 ' 



( - 1 
1 -1 



(3.1) 



(3.2) 



B 



2Ax 



1 -1 
1 / 



and solution 



The exact amplification factor (9 — kAx) 



V-i 

u{t + At) = e Am u{t). 



(3.3) 



-%t] sin e 



with 



jAt 



(3.4) 



(3.5) 



is unitary and causes a phase-shift of each Fourier component. In the limit of Ax — > 0, the phase-shift becomes uniform 
for all the Fourier components e lkx — » Q lk i x - vAt ) ^ resulting in a uniform shift of the entire function u(x) — > u(x — vAt), 
which is the exact solution to (|3.ip . Any Taylor expansion of (|3.4[) will produce algorithms with a non- unitary g, 
resulting in unwanted dissipations or instability. The situation here is much more delicate than in the diffusion case. 
The natural decomposition is similarly, 



iV 



(3.6) 



where 



B, 



2Ax 



It follows that 



-1 

1 



with B 



AT 



2Ax 



( o 



1\ 



0/ 



(3.7) 



c — s 
s c 



,AtB r 



(3.8) 
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where now 

c = cos(r//2) and s = sin(??/2). (3.9) 

Each e AtBj only updates Uj and Uj + \ as 



Uj = cuj — suj + i 
u' j+1 = suj + cuj+i. (3.10) 

As in the diffusion case, the above updating can be recasted into the following forms for algorithms 1A and IB, with 
1A given by 

u\ = CU\ — SU 2 

u' 2 = SU* + 1*2 — Sl*3 

u'j = su' J _ 1 + uj - su 3+1 (2 < j < N) (3.11) 
u' N — su' N _ 1 + u N — su* 



cu\ = su' N + cu\ — SU2- 



and IB given by 



u\ — sun + cui 

u' N = SMJV-1 + Upf — su* 

u'j = suj-! + Uj - su' J+1 (2<j< N) (3.12) 
u' 2 = su* + U2 — su' 3 
cu\ = sun + cu\ — su 2 . 

In contrast to the diffusion case, these algorithms are not exactly norm-preserving for periodic boundary condition. 
By adding up both sides of the above algorithms, one finds that what is preserved by 1A is not the usual norm 
N = Y^j=i u i-> but a modified norm given by 

Nia = N + (-^ - l) Ul (3.13) 



Similarly, what is preserved by IB is 



N 1B = N + (— (3.14) 

1 + s 



If initially u\ = 0, then TVia = Nib — No, where No is the initial norm. As the system evolves, each algorithm's 
actual norm will evolve as 

N 1A = N - (-^— - l)ui, 
1 — s 

N 1B = Nq — {t—7 (3.15) 

1 + s 

The error is due to a single point u\, where it is the only point not updated twice immediately. As the wave form 
travels around the periodic box, u\ will trace out the shape of the wave and imprint that as the error of the norm 
in time. For a sharp pulse, the norm error will return to zero after the pulse peak has passed through u\. Thus 
norm-preservation will be periodic. For the advcction equation, this is a small effect, and is secondary to the phase 
and oscillation error mentioned below. However, this error will be important in the next section. 

If the boundary values u[, u' 2 and u' N are ignored for now, then again the resulting second-order algorithm is unique, 
independent of the order of applying 1A or IB. The amplification factors are all unitary: 

1 — se i0 

9ia = z — 7) = cxp(-i0i A ), (3.16) 

1 — se ia 



3 1B = 1 r—fe =exp(-i(ftiB,) (3.17) 



1 + se 
1 + se 

g 2 = g 1B (At/2)g 1A (At/2) 
1 - i2{s/c 2 ) sin 



1 + i2(s7S 2 )sin6> 



cx P H</> 2 ), (3.18) 
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with phase angles 

_-, / s sin 9 

4>ia = 2 tan 



'IB 



2 tan" 



1 — s cos ( 
if s sin 9 



1 + s cos ( 
02 = 0ia(^/2) + 0ib(^/2), 

= 2 tan" 1 ^— Hf^sin^ , (3.19) 

where here 

s = sin(r//4) and c = cos(r?/4). (3.20) 

Since (Jia and <?ib are not complex conjugate of each other, their phase errors do not exactly cancel. Their residual 
difference is the error of the second-order algorithm. 

Algorithms (|3.11l) and (|3.12l) are the corresponding Saul'yev's schemes for solving the advection equation. The 
coefficient here is s — sin(r//2) rather than Saul'yev's coefficient of s = rj/2. This explains why it makes no sense to 
apply Saul'yev's schemes at s > 1, since they can no longer be derived from the fundamental updating matrix (|3.10[) 
with a real c = vl — s 2 . At s > 1, Saul'yev's schemes are in fact unstable, suffering from spatial amplification^, 
despite the unimodulus appearance of (|3.16l) and (|3.17p . This is easy to see in the case of algorithm 1A. If initially 
Uj — for j > J, but uj_i 7^ 0, then according to (|3.11|1 . u' J+n = s n+1 u' J _ l increases without bound as a function of 
n. Even the case of s = 1 is pathological. For Saul'yev's coefficient s = r\/2 = 1, one has 

91A = -o lkAx = -e ik ^ At and g 1B = e~ lkAx = e~^ v ^ At . (3.21) 

Under algorithm 1A, Fourier mode e lkx will flip its sign and propagate with velocity —v/2. Under IB, it will propagate 
with velocity v/2. The resulting second order algorithm then leaves the Fourier mode stationary with only a sign flip. 
This is completely contrary to the behavior of the exact solution and is a source of great error for Saul'yev's schemes. 
As we will show below, alternative choices for s will eliminate such unphysical behaviors. 

While the derived choice of s — sin(r)/2) is unconditionally stable for all ry, the resulting algorithms 1A and IB have 
huge phase errors, and are no better than Saul'yev's choice of s = rj/2. This is because in comparison with the exact 
phase angle, 

cj) ex =7 ] sm9 = r ] 9- ^9 3 + ■ ■ ■ (3.22) 
6 



algorithms 1A and IB have expansions 



23 8(1 + 8) 3 

2s s(l - s) 3 

6iB - i+- s 9 -3(rTsr e + 



(3.23) 



and neither s = rj/2 nor s — sin(ry/2) can result in a first-order coefficient of 9 matching that of 4> ex exactly. The 
choices of s that can do this are, for 1A, 

2S V -> s = -^-, (3.24) 



1 - s 2 + 7/ 

and for IB, 

2s 



V s=-^. (3.25) 

1 + s 2 — i] 

This then reproduces the Roberts and Weis a 23 > 24 forms of the Saul'yev-type algorithm and will be denoted as RW1A 
and RW1B. For r) > 0, only RW1A is unconditionally stable and RW1B is limited by spatial amplification to rj < 1. 
The pathological behavior of 1A at s = 1 can no longer occur at any finite 77. 
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For the above choices of s, the corresponding phase angles are 

fcA = ^-(6 + T + 12 } 



= v0-(l~^- + ^)e 3 + ---, (3.26) 



and the modified norms are 



N 1A = N + {y/T+7j - l)tii, (3.27) 
JVib = iV+( v /T _ ^-l)tti. (3.28) 



The second order algorithm from concatenating RW1A and RW1B is 

4>2 = &a(»?/2) + &b(»7/2) 



= ^ _ (1 + ??)e 3 + (JL + + JlL)qs + . . . (3 29) 

' V 6 48' v 120 192 1280 ; V ' 

This second-order advection algorithm will be denoted as RW2. Because RW1B is limited by spatial amplification to 
rj < 1, RW2 is limited in stability to rj < 2. 

To generate a stable second-order algorithm for all r/, one can concatenate 1A and IB with the same s. To match 
g ex to first order in 6 then requires 

rh + rh^' - tS = I - '-li/i+^-'l- (««) 

The resulting amplification factor is, according to (|3.18|) . 

_ 1-1(77/2) sing 

92 " 1 + ^/2)81110' (3 ' 31) 




which is precisely the implicit Crank-Nicolson amplification factor. Since by (I3.30[) , s < 1, and c = vl — s 2 is well- 
defined for all 77, the algorithm is unconditionally stable and can be applied to periodic boundary problems via the 
the fundamental updating (13.101) . Corresponding to (|3.31[) . the phase-angle has the characteristic expansion, 

^ 2 = r} e-( r l + 1 f-)e s + (^L + r t + ^t)e b + --- (3.32) 

V2 1 v 6 12^ 420 24 240 ' V ' 

= 7 1 sm8 + i 1 3 F 3 (6)+ V 5 F 5 {6) + --- (3.33) 

where now the time parameter is 77 and the original "Hamiltonian" is h — sin 8. We shall designate this second-order 
algorithm, with s given by (|3.30p . as A2C. The second-order algorithm corresponding to Saul'pev's choice of s = rj/4 
will be denoted as A2S, and the initially derived result of s — sin(?7/4) as A2. 

On the left of FigJU the phase error of these symplectic algorithms are exaggerated and compared to the explicit 
but dissipative Lax-Wendroff (LW) scheme at a large value of 77 = 0.7. The original 1A and IB algorithms have huge 
phase errors but are mostly cancelled in the second-order algorithm A2. Even so, algorithm A2's error curve has a 
finite slope at 9 — 0, as shown on the right of Fig0J By construction, schemes RW1A, RW1B, RW2, A2C and LW all 
have zero error slopes at = 0. This is a crucial advantage of A2C over A2. Also, A2C is stable for all 77, while RW2 
is limited by spatial amplification to 77 < 2. 

The importance of having a zero error slope in the phase angle can be better appreciated from the following 
considerations. Let 

, . f xuix, t)dx 
X(t) = \ / ■ (3.34) 

J u{x, t)dx 

Since the solution to the advection equation is u(x, t) — uq{x — vt), we have the exact result 

t A .x J xu Q (x - vAt)dx J(y + vAt)u (y)dy 

x(At) = — r — - ( T-rr- = f — l~T~j = X(0) + vAt. (3.35) 

J uo{x — vAtfdx J uo(yjdy 
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FIG. 4: The phase error of various advection algorithms at r\ = 0.7. Left: The unmodified symplectic algorithms are denoted 
as 1A, IB, and A2. The Robert- Weiss versions are denoted as RW1A, RWIB and RW2. A2C is the second-order algorithm 
with the Crank-Nicolson amplification factor. LW is the Lax-Wendroff scheme included for comparison. Right: The phase 
errors of fourth and sixth-order algorithms composed out of three, five and seven second-order algorithms Ti. Their phase 
errors are compared to that of running the second-order algorithm three, five and seven times at reduced time steps. The T2 
used here is A2C. Previous second-order algorithms are also included for a close-up comparison. 



Multiplying algorithm 1A (|3.1ip by j and sum over j yields 



N N N N 

N N N N 

j=i i=i i=i j=i 

For a localized pulse far from the boundary, the norm can be consider conserved, 

N N 

One then has the discrete version of (|3.35l) 



"3 _ i-ij- 



which will reproduce the displacement exactly if 



2^7=1 



1 - s 



■Ax 



2s vAt 



1 - s Ax 



(3.36) 



(3.37) 



(3.38) 



(3.39) 



which is the condition (|3.24p for a zero error-slope. Similarly for IB satisfying (|3. 251) . Far from the boundary, 
symplectic advection algorithms with a zero-error slope in the phase angle would exactly preserve the first two 
moments of (x n ). 

In FigjS] we show the working of these algorithms in propagating an initial profile 



u(x, 0) = exp 



(§)' 



(3.40) 
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The power of 6 was chosen to provide a steep, but continuous profile so that both the phase error and the oscillation 
error are visible. If the profile were too steep, like that of a square wave, the oscillation error would have overwhelmed 
the calculation before the phase error can be seen. The oscillation errors in all these symplectic algorithms are 
primarily due to the oscillation error in algorithm 1A. Algorithm IB has a much smaller oscillation error. Because all 
the algorithms are essentially norm-preserving, oscillation errors is inherent to any scheme which does not preserve 
the positivity of the solution^!. 




FIG. 5: The propagation of initial profile (|3.40p fives times around a periodic box of [-10,10] with Ax = 0.025, At = 0.02, 
v = 1, and r\ = 0.8, corresponding to 5000 iterations of each algorithm. If there were no phase error, the profile would remain 
centered on x = 0. Algorithms A2 and A2S have large and positive, phase errors. The oscillation errors in these second-order 
symplectic algorithms are predominately due to the imbedded 1A algorithm. Algorithm RW1B has a much smaller oscillation 
error and is comparable to the oscillation error in the dissipative Lax-Wendroff scheme (bright green line) . 

Higher order advection algorithms can again be constructed by the method of composition (|1.17[) . Since g 2 is 
unitary, any product of g 2 is also unitary. Thus to preserve unitarity, one must use only a single product composition, 
rather than a multi-product expansion as in the diffusion case. (However, as noted in the last secrion, this violation 
of unitarity in MPE is small with increasing order. At sufficiently high order, this violation is beyond machine 
precision and is indistinguishable from a truly unitary algorithm—. For simplicity, we will only consider strictly 
unitary algorithms in this discussion.) For a single product composition, the resulting phase angle is just a sum of 
</> 2 's. The simplest fourth-order composition, the Forest- Ruth (FR) algorithm^— is given by 

T[ R (At) = T 2 (a 1 Z\i)T 2 (a ^)T 2 (a 1 Z\i), (3.41) 

with ai = 1/(2 — b), ao = —b/{2 — b) and b = 2 1 / 3 . The coefficients ao and a\ satisfy the consistency condition 
2ai + ao ~ 1 and the fourth-order condition 2a\ + Oq = 0. If we take T 2 (Z\i) to be A2C, then the phase angle for 
T[ R (At) is just (|3~32j) with all rf terms removed, 

0£R = rfi - V + ( -5- - 5.29145^^) 9 5 + ■ ■ ■ , (3.42) 
V4 1 6 V 120 240/ ' V ' 

which is then correct to fourth-order in 9. This is not true if we take T 2 (Z\t) to be the original algorithm A2. That 
fourth-order time- marching algorithm's phase angle will still have a small error slope at 9 — 0. One should therefore 
only uses A2C to compose higher order algorithms. 
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The large numerical coefficient in (I3.42[) is due to 2a\ + Oq = —5.29145, reflecting the fact that FR has a rather 
large residual error. A better fourth-order algorithm advocated by Suzuki^ (S4) at the expense of two more T 2 is 

Tf = T 2 (a 1 At)T 2 (a 1 At)T 2 (aoAt)T 2 {a 1 At)T 2 (a 1 At) (3.43) 

where now ai — 1/(4 — 4 1 / 3 ), clq = — 4 1 / 3 a 1 , and Aa\ + 0Sq = —0.074376, which is nearly sixth-order. 
At the expense of two more T 2 , one can achieve sixth-order via Yoshida's algorithm^ (Y6), 

Tj = T 2 ( a3 Z\<)T 2 (a 2 Z\i)r 2 (a 1 Z\i)T 2 ( ao Z\<)T 2 (a 1 Z\i)7 1 2(a 2 ^)? 1 2(fl3^) (3.44) 

with coefficients 

ai = -1.17767998417887 a 2 = 0.235573213359357 

a 3 = 0.784513610477560 and a = 1 - 2(d + c 2 + c 3 ). (3.45) 

For eighth and higher order algorithms, see RefsUH andl25l 

The phase errors of these higher order algorithms at 77 = 0.7 are shown at the right of Fig|4] Since each algorithm 
applies T 2 (taken to be A2C) n (=3,5,7) times, they are compared to the phase error of running T 2 n times at a 
reduce time step of At/n. Algorithms S4 and Y6 beat their target comparisons by orders of magnitude. There is a 
clear advantage in going to higher-order algorithms for solving the advection equation. 




FIG. 6: Comparing the convergence of various unconditionally stable symplectic advection algorithms. The solid lines are power 
laws of the form a + bAt n seeking to verify the order of the algorithm. Higher order algorithms composed of n second-order 
algorithms A2C are compared to n*A2C at a reduce time-step size of At/n. FR, S4 and Y6 requires 3, 5 and 7 runs of A2C 
respectively. In this computation of (|3.46[) . all algorithms met or exceeded their nominal order of convergence. See text for 
details. 

In Fig|6l the convergence of these higher order algorithms are compared. The range of the time steps used, 
At = 0.01 — 0.20, corresponds to 77 = 0.4 — 8.0. The same profile (|3.40[) is initially centered at x = —5 and propagated 
to x = 5 at v = 1. The time steps are chosen as At = 10/m, so that m iterations exactly give t = 10. Since all 
algorithms satisfy (|3.35[) despite the oscillation errors, we compute the expectation value 

SjliCWM ... 

({x}) = 3 N (3.46) 
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with respect to the absolute value of the propagated profile. The phase errors for A2 and A2S are known to be large 
from Fig|5J and are not included in this comparison. The solid lines are fitted power laws of the form a + bAt n . The 
first suprise is that RWlA's result cannot be fitted with n = 1. The fitted line is a fit with n = 3/2. The algorithm 
A2C can be well fitted with n = 2. This is specially clear in the case where A2C is applied seven times at step 
size At/7. The fourth-order Forest-Ruth (FR) and Suzuki (S4) algorithms can only be fitted with n = 5, and the 
sixth-order Yoshida (Y6) algorithm with n = 7. They all converged to a value of a = 4.999438 which is below the 
exact value of 5. This is related to the grid size error. Halving the grid size to Ax — 0.0125 gives a = 4.999997. 



IV. SYMPLECTIC ADVECTION-DIFFUSION ALGORITHMS 



The advection-diffusion equation, 



has the exact operator solution 



du du jj^ 2u 

dt dx dx 2 ' 



u{x,At)=e- vAt £ +DM -&u(x,0). 



(4.1) 



(4.2) 



If v and D are just constants, then since 



r_9_ 

L dx 1 dx z 



] = 0, one has 



u{x,At) 



-vAt-g- 



•u(x, 0) 



u{x, At) 
u(x - vAt,Ai), 



(4.3) 



where u(x, At) is the diffused solution. The complete solution is therefore the exact diffused solution u(x, At) displaced 
by vAt. 

For periodic boundary condition, our matrices also commute, [A, B] = 0, so that the discretized version also holds, 



u(t + At) = e AtA e AtB u(t) 



(4.4) 



Thus arbitrary high order algorithms can be obtained by applying higher order advection and diffusion algorithms in 
turns from the previous sections. In the case of spatially dependent D(x) or v(x) where [A, B] ^ 0, one can do the 
second-order splitting 



u{t + At) = e^ tA e 4 * B e^ tA u(i) 



(4.5) 



and apply higher order MPE algorithms. However, for periodic boundary condition, this way of solving the advection- 
diffusion equation cannot conserve the norm. Consider the case of applying the advection algorithms RW1A, RW1B 
followed by any norm-conserving diffusion algorithm. From (|3.27p . the change in the modified norm after At would 
be 



Nia = (V^+V - 1)K - ui). 



(4.6) 



For RW1A, u\ is a discontinuous point higher than its adjacent neighbors u n and u 2 - Consequently, after the diffusion 
step, Ui < u\ and there is a loss of normalization in (I4.6[) . For RW1B, u\ is a discontinuous point lower than its 
adjacent neighbors u n and ui- After the diffusion step, u\ > u\, and again results in a loss of normalization: 



n' ib - N m = (VI - ^ - i)K - m). 



(4.7) 



As will be shown, this loss is small for small D, but is irreversible and accumulative after each orbit around the 
periodic box. For fixed boundary with u\ = 0, there is no such norm-conserving problem. 

An alternative is to update the advection and diffuion steps simultaneously. In this case, one might decomposing 
A + B into a sum of 2 x 2 matrices as done previously, 



A, 



B, 



D 

A? 



-1 1 
1 -1 



2Ax 



\ 



-1 

1 



(4.8) 
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resulting in 

e AtC > = | " ^ I , e Atc » = - n | , (4.9) 

with 

a = e - r cosh^, /3 = e- r (r + V /2)^^, X = e~ r (r - r]/2)^^-^-, (4.10) 

ip ip 





and ip = \Jr 2 — (r//2) 2 . The corresponding Saul'pev form of the 1A algorithm is then 

u'j = /3u' j _ 1 + juj + Xuj + i (4.11) 

where 7 = a 2 — f3X = e~ 2r remains the determinant of the updating matrix. However, for the Saul'pev form (|4.11l) to 
be norm-preserving, one must have 

/3 + 7 + A = l. (4.12) 

Surprisingly, this is grossly violated by (|4.10|) when both r and 77 are non- vanishing. 

As we have learned in the previous two sections, any such initial algorithm can be far from optimal. Therefore, one 
may as well begin with an assumed updating matrix, 

u'j = au,j + Xuj+i 
u'j +1 — (3uj + auj + i (4-13) 

and determine its elements by enforcing norm-conserving condition (|4.12l) and by matching the expansion coefficients 
of the exact amplification factor. The sequential applications of this updating matrix yields Saul'pev-type algorithms 
1A (|4~TTJ) and IB, 

u'j = fiuj-i + jUj + Xu' j+1 (4.14) 



The determinant 7 = a 2 — (3X is to be regarded as fixing a as a function of 7 and /3 via a = y/j + PX. The resulting 
amplification factors are then 

The norm condition (j4.12[) fixes A in terms of /3 and 7. In terms of 7 and /3 algorithms 1A and IB have expansions, 

hm . g^l ^y ^Oyi (4-16) 

j + /3 2(7 + py 

Matching the first and second order coefficients of the exact exponent 

h ex = ir]sm(8) + 4rsin(6»/2) 2 

= ir/9 + r9 2 + O{0 3 ), (4.17) 

then completely determines, for 1A and IB respectively, 

I-7 + 17 1- wr 2 

2 + 77 1 + wr 2 + 77(3 + 77) 

1 — 7 + 777 1 — wr 2 

/?= 9 r r/ 7=T-— w = - R y 4.19 

2 — 7/ 1 + wr 2 — 77(3 — 17) 
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These are the generalized Roberts- Weiss algorithms for the advection-diffusion equation. 
For any choice of /? and 7, the modified norm including the boundary effect now reads 



N 1A = N+( 
N 1B = N + ( 



1-/3 



l)t*i 



1-A 

Remarkably, for the above generalized RW algorithms, one has 

and 



-l)ui 



1-/3 



V 



1-A 



(4.20) 
(4.21) 

(4.22) 



The modified norms (I4.20[) and (|4.21[) are therefore the same as the pure advection cases of (|3.27[) and (I3.28[) . For 
these two first order advection-diffusion algorithms, their norm-conservation in a periodic box will then be periodic, 
as in the pure advection case. This is shown in FigJTl However, as soon as one concatenate them into algorithm RW2, 
the loss of norm is irreversible. This is because each algorithm will behave as a diffusion algorithm for the other. The 
norm-loss mechanism described earlier will then apply. 



0.01 



0.00 



-0.01 - 



-0.02 




FIG. 7: The normalization error of various algorithms when propagating a Guassian profile in a periodic box of [0,10] with 
Ax = 0.05, At = 0.033, v = 1, D = 0.005, r\ = 0.66 and r = 0.066. The profile is initially centered at x = 5. At t = 5, 15, 25, 25, 
the Gaussian peak is at the edge of the periodic box. First-order advection-diffusion algorithms RW1A and RW1B conserve the 
norm periodically. All higher than first-order algorithms suffer loss of normalization irreversibly, though very small for MPE 
algorithm T4 and A/D. The latter is applying the advection algorithm A2C and the diffusion algorithm D2S sequentially. 



The second-order algorithm's amplification factor is 

/ 7 + Ae 

92 ~ 




7 + Pe~ i6 



v 1 - fie-* 9 I \ 1 - Ae+ 
where 7 = r y(At/2), etc.. In terms of 7 and (3, hi has the expansion, 

/ (I + 7X7- 1 + 2^ , , 



0(9'). 



(4.23) 



(4.24) 
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Matching this to the first order coefficient of the exact exponent (|4. 1 7[) determines 



and 



^(1-7) + ^(1 + 7)3; (4-25) 



A=i(l-7)-i(l + 7)s. (4.26) 



where s has been previously defined by p. 301) . In terms of only 7 

, (1-7) (1 + 3S 2 
'(1+7) 

and matching the second order coefficient in (j4.17[) determines 



W + 2 Ti , ~( ^ ~n 2 g2 + °( g3 )> ( 4 - 27 ) 



^ = 1W2 w . th u; = (1 _ s2) 2 /(1 + 3s2) _ (428) 
1 + ior/2 

If r\ = 0, s^= 0, one recovers (|2.33l) . which is the second-order diffusion algorithm D2S. If r = 0, then 7 = 1 and one 
recovers the second-order advection algorithm A2C with A = —/3 = —s. We shall refer to this second-order algorithm 
as AD2C. AD2C is an unconditionally stable algorithm which requires only half the effort of algorithm A/D, which 
applies A2C and D2S scquentually. However, AD2C has a greater irreversible norm-error when applied to periodic 
boundary problems. This is shown in Fig|7] 




FIG. 8: The propagation of a Guassian profile in a periodic box of [0,10] with Ax = 0.05, At = 0.033, v = 1, D = 0.1, r\ = 0.66 
and r = 1.33. The profile is initially centered at x — 5. All three profiles produced by algorithms RW1A, AD2C and A/D are 
in essential agreement prior to the pulse peak hitting the right periodic edge. As the profiles reappear from the left, the norm 
of AD2C is noticeably lower. 

This norm-error for periodic boundary condition can be greatly reduced by going to higher orders. FiglT] shows the 
result for the fourth-order MPE algorithm T 4 , with AD2C as T 2 . For r small, surprisingly, even the negative-coefficient 
algorithm FR is stable, but with an error comparable to second-order algorithms. 
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In Fig J8] we illustrate the effect of this norm-loss error at a large value of r = 1.33. For clarity, only results from 
three representative algorithms are shown. Algorithms A/D and RW1A have small or only periodic norm-losses and 
remained in agreement after the Gaussian peak has reappeared from the left. However, algorithm AD2C suffers an 
irreversible norm-loss and its peak is noticeably lower. 



V. CONCLUDING SUMMARY 

In this work, we have shown that explicit symplectic finite-difference methods can be derived in the same way as 
symplectic integrators by exponential splittings. The resulting sequential updating algorithms reproduce Saul'pev's 
unconditionally stable schemes, but is more general and can be applied to periodic boundary problems. In contrast 
to Saul'pev's original approach, where the algorithm is fixed by its derivation, symplectic algorithms can be system- 
atically improved by matching the algorithm's amplification factor more closely to the amplification factor of the 
semi-discrctizcd equation. One key contribution of this work is the recognition that, for finite difference schemes, 
their amplification factors should be compared, not to the continuum growth factor, but to the amplification factor 
of the semi-discretized equation. The exponent of this amplification factor then serve as the "Hamiltonian" for de- 
veloping symplectic finite-difference algorithms. By requiring the algorithm's modified "Hamiltonian" to match the 
original "Hamiltonian" to the leading order, one produces all known, non-pathological first-order Saul'pev schemes 
and many new second-order algorithms for solving the diffusion and the advection equation. As a consequence of this 
formal correspondence with symplectic integrators, existing methods of generating higher order integrators can be 
immediately used to produce higher order finite-difference schemes. 

The generalization to higher dimensions can be done by dimensional splitting, resulting in unconditionally stable, 
alternate-direction-explicit methods. The generalization to non-constant diffusion and advection coefficients is a topic 
suitable for a future study. The coefficients must frozen in such a way that one can recover Saul'pev's asymmetric 
schemes from their more basic sequentual updatings. 

Acknowledgments 

This work is supported in part by the Austrian FWF grant P21924 and the Qatar National Research Fund (QNRF) 
National Priority Research Project (NPRP) grant # 5-674-1-114. I thank my colleague Eckhard Krotscheck and the 
Institute for Theoretical Physics at the Johannes Kepler Univeristy, Linz, Austria, for their wonderful hospitality 
during the summers of 2010-2012. 



1 R. Courant, K.O. Friedrichs, H. Lewy "Uber die partiellen Differenzengleichungen der mathematischen Physik", Maht. Anal. 
100, 32-74 (1928). 

2 V. K. Saul'yev, "On a method of numerical integration of a diffusion equation", Dokl. Akad. Nauk. SSSR, (in Russian) 115, 
1077-1079 (1957). 

3 V. K. Saul'yev, Integration of equation of parabolic type by the nethod of nets, Pergamon Press, New York, 1964. 

4 B. K. Larkin, "Some stabel explicit difference approximation to the diffusion equation", Math. Comput. 18, 196-201 (1964). 

5 D. J. Evans and A. R. B. Abdullah, "Group Explicit Methods for Parabolic Equations", Intern. J. Computer Math. 14, 
73-105 (1983). 

6 D. J. Evans, "Alternating group explicit methods the diffusion equations", App. Math. Modelling, 9, 201-206 (1985). 

7 M. Creutz and A. Gocksch, "Higher-order hydrid Monte-Carlo algorithms", Phys. Rev. Letts. 63, 9 (1989). 

8 E. Forest and R. D. Ruth, "4th-order symplectic integration", Physica D 43, 105 (1990). 

9 H. Yoshida, "Construction of higher order symplectic integrators", Phys. Lett. A150, 262-268, (1990). 

10 M. Suzuki, "Hybrid exponential product formulas for unbounded operators with possible applications to Monte Carlo simu- 
lations", Phys. Lett. A 146, 319 (1990). 

11 H. Yoshida, "Recent progress in the theory and application of symplectic integrators", Celest. Mech. Dyn. Astron. 56, 27 
(1993). 

12 E. Hairer, C. Lubich, and G. Wanner, Geometric Numerical Integration, Springer- Verlag, Berlin-New York, 2002. 

13 R. I. McLachlan and G. R. W. Quispel, "Splitting methods", Acta Numerica 11, 241 (2002). 

14 H. F. Trotter, "Approximation of semi-groups of operators", Pacific J. Math. 8, 887-919 (1958) 

G. Strang, "On the construction and comparison of difference schemes" SIAM J. Numer. Anal. 5, 506-517 (1968). 

16 Q. Sheng, "Solving linear partial differential equations by exponential splitting", IMA Journal of numberical anaysis, 9, 
199-212 (1989). 

17 S. A. Chin, "Multi-product splitting and Runge-Kutta-Nystrom integrators", Cele. Mech. Dyn. Astron. 106, 391-406 (2010). 



22 



S. A. Chin and Jurgen Geiser, "Multi-product operator splitting as a general method of solving autonomous and nonau- 

tonomous equations", IMA Journal of Numerical Analysis, 31 1552-1577 (2011); doi: 10.1093/imanum/drq022 

M. Schatzman, "Numerical integration of reaction-diffusion systems", Numerical Algorithms 31 247-269 (2002). 

S. Blanes, F. Casas and J. Ros, Extrapolation of symplectic integrators, Celest. Mech. Dyn. Astron., 75 (1999)149-161 

W. Hundsdorfer and J. G. Verwer, Numerical Solution of Time- Dependent Advection-Diffusion-Reaction Equations, P.119, 

Springer- Verlag Berlin Heidelberg 2003. 

R. E. Zillich, J. M. Mayrhofer and S. A. Chin, "Extrapolated high-order propagators for path integral Monte Carlo simula- 
tions", J. Chem. Phys. 132, 044103 (2010). 

K. V. Robert and N. O. Weiss, "Convective difference schemes", Math. Comput. 20, 272-299 (1966) 

L. J. Campbell and B. Yin, "On the stability of Alternating-Direction Explicit Methods for Advection-Diffusion Equations", 
Num. Methods for PDE 23, 1429-1444 (2007); DOI: 10.1002/num.20233 

R. I. McLachlan, "On the numerical integration of ordinary differential equations by symmetric composition methods", 
SIAM J. Sci. Comput. 16, 151 (1995). 

S. Blanes and P. C. Moan, "Practical symplectic partition Runge-Kutta menthods and Runge-Kutta Nystrom methods", J. 
Comput. Appl. Math. 142, 313 (2002). 



